#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Merge & Align (Iterative)                                          #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 01/03/2007                                             #
# Last Modification: 01/03/2007                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 40
#
# MANUAL: This script allows refining the data of the selected images.
#
# MANUAL: All selected images are merged together, and subsequently aligned against the merged average. This can be done iteratively a preset number of times. Double-click onto the iteration counter after the script name, to adjust that number, if needed.
#
# PUBLICATION: 3D reconstruction of two-dimensional crystals: <A HREF="http://www.ncbi.nlm.nih.gov/pubmed/26093179">Arch Biochem Biophys 581, 68-77 (2015)</A>
# PUBLICATION: 3D Reconstruction from 2D crystal image and diffraction data: <A HREF="http://dx.doi.org/10.1016/S0076-6879(10)82004-X">Methods Enzymol. 482, Chapter 4, 101-129 (2010)</A>
# PUBLICATION: 2dx - Automated 3D structure reconstruction from 2D crystal data: <A HREF="http://journals.cambridge.org/action/displayAbstract?aid=1943200">Microscopy and Microanalysis 14(Suppl. 2), 1290-1291 (2008)</A>
# PUBLICATION: 2dx_merge - Data management and merging for 2D crystal images: <A HREF="http://dx.doi.org/10.1016/j.jsb.2007.09.011">J. Struct. Biol. 160(3), 375-384 (2007)</A>
#
# DISPLAY: MergeResolution
# DISPLAY: zstarwin
# DISPLAY: zstarrange_real
# DISPLAY: RESMIN
# DISPLAY: RESMAX
# DISPLAY: merge_res_limit
# DISPLAY: tempkeep
# DISPLAY: realcell
# DISPLAY: realang
# DISPLAY: ALAT
# DISPLAY: MergeStepSize
# DISPLAY: IBOXPHS
# DISPLAY: SYM
# DISPLAY: avrgamphsNUMBER
# DISPLAY: avrgamphsRESOL
# DISPLAY: MergeIQMAX
# DISPLAY: Merge_Reference_IQMAX
# DISPLAY: MergeHKMAX
# DISPLAY: Merge_Reference_HKMAX
# DISPLAY: AMP_Scale_Factor
# DISPLAY: ILIST
# DISPLAY: refbeamtilt
# DISPLAY: reftiltgeo
# DISPLAY: merge_reference
# DISPLAY: merge_ref_num
# DISPLAY: merge_comment_1
# DISPLAY: merge_comment_2
# DISPLAY: merge_comment_3
# DISPLAY: merge_comment_4
# DISPLAY: merge_comment_5
# DISPLAY: merge_comment_6
# DISPLAY: merge_comment_7
# DISPLAY: merge_comment_8
# DISPLAY: merge_comment_9
# DISPLAY: merge_refine_iterations
# DISPLAY: max_amp_correction
# DISPLAY: make_reference
# DISPLAY: merge_data_type
# DISPLAY: plotres_rings
# DISPLAY: scalimamp3d_rref
# DISPLAY: scalimamp3d_BXYMINMAX
# DISPLAY: scalimamp3d_BZMINMAX
# DISPLAY: scalimamp3d_BEXTRA
# DISPLAY: resolutionplot_RESMAX
# DISPLAY: resolutionplot_bins
# DISPLAY: RFACAMP
# DISPLAY: Thread_Number
# DISPLAY: RESMIN
# DISPLAY: RESMAX
# DISPLAY: merge_alsoevenodd
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
#
set ALAT = ""
set AMP_Scale_Factor = ""
set avrgamphsNUMBER = ""
set avrgamphsRESOL = ""
set Calc_from_zstarrange = ""
set IBOXPHS = ""
set ILIST = ""
set lattice = ""
set max_amp_correction = ""
set merge_data_type = ""
set merge_ref_num = ""
set merge_reference = ""
set Merge_Reference_HKMAX = ""
set Merge_Reference_IQMAX = ""
set merge_refine_iterations = ""
set merge_res_limit = ""
set MergeHKMAX = ""
set MergeIQMAX = ""
set MergeResolution = ""
set MergeRminRmax = ""
set MergeStepSize = ""
set num_reflections_fitted = ""
set num_reflections_FOM50 = ""
set plotres_rings = ""
set realang = ""
set realcell = ""
set refbeamtilt = ""
set reftiltgeo = ""
set RESMAX = ""
set RESMIN = ""
set resolutionplot_bins = ""
set resolutionplot_RESMAX = ""
set RFACAMP = ""
set scalimamp3d_BEXTRA = ""
set scalimamp3d_BXYMINMAX = ""
set scalimamp3d_BZMINMAX = ""
set scalimamp3d_rref = ""
set SYM = ""
set tempkeep = ""
set Thread_Number = ""
set zstarrange = ""
set zstarrange_real = ""
set zstarwin = ""
set ctfrev = ""
set calculate_subvolume = ""
set number_of_beads = ""
set number_refinement_iterations = ""
set membrane_height = ""
set maximum_amplitude_refinement = ""
set merge_alsoevenodd = ""
#
#$end_vars
#
setenv OMP_NUM_THREADS ${Thread_Number}
#
set ccp4_setup = 'y'
source ${proc_2dx}/initialize
#
set ITAXASTEP = "1"
set RTAXASIZE = "2.0"
set ITANGLSTEP = "1"
set RTANGLSIZE = "0.5"
#
set split = ($realcell:as/,/ /)
set cellx = $split[1]
set celly = $split[2]
#
echo "cellx = ${cellx}"
echo "celly = ${celly}"
echo "cellz = ${ALAT}"
#
set cellxm1 = `echo ${cellx} | awk '{ s = $1 - 1 } END {print s}'`
set cellym1 = `echo ${celly} | awk '{ s = $1 - 1 } END {print s}'`
set ALATm1 = `echo ${ALAT} | awk '{ s = $1 - 1 } END {print s}'`
#
set dirfile = "2dx_merge_dirfile.dat"
set dirnum = `cat ${dirfile} | wc -l`
set maxthread = `echo ${Thread_Number} ${dirnum} | awk '{if ($1<$2/2) { s = $1 } else { s = int($2 / 2) }} END { print s }'`
if ( ${maxthread} == "0" ) then
  set maxthread = 1
endif
#
set number = 1
if ( ${ILIST} == "n" ) then
  set IVERBOSE = 1
else
  set IVERBOSE = 6
endif
#
set scriptname = 2dx_mergeRefine
set merge_modus="3D"
#
\rm -f LOGS/${scriptname}.results
#
# The following is to make sure that for the next "Import Images", the default is correctly initialized.
set initialization_reset = "y"
set initialization_executable = "y"
echo "set initialization_reset = ${initialization_reset}" >> LOGS/${scriptname}.results
echo "set initialization_executable = ${initialization_executable}" >> LOGS/${scriptname}.results
#
if ( ${merge_ref_num} != "0" ) then
  ${proc_2dx}/linblock "ERROR: Iterative refinement does not make sense with"
  ${proc_2dx}/linblock "a saved static reference."
  ${proc_2dx}/linblock "Use last merge result as reference (register 0),"
  ${proc_2dx}/linblock "or use the REFINE script instead."
  ${proc_2dx}/protest "Aborting."
endif
#
set zstarrange = `echo ${zstarrange_real} ${ALAT} | awk '{ s = 1.0 / ( $1 ) } END { print s }'`
echo "set zstarrange = ${zstarrange}" >> LOGS/${scriptname}.results
${proc_2dx}/linblock "Calculating zstarrange as ${zstarrange} reciprocal Angstroems."
#
set zmin = `echo ${zstarrange} | awk '{s = -$1} END {print s}'`
set zminmax = `echo ${zmin},${zstarrange}`
echo zminmax = ${zminmax}
# 
set zstarwintmp = `echo ${ALAT} | awk '{ s = 1.0 / ( 2.0 * $1 ) } END { print s }'`
set zstarwinbad = `echo ${zstarwin} ${zstarwintmp} | awk '{ if ( abs ( $1 - $2 ) > 0.001 ) { s = 1 } else { s = 0 } } END { print s }'`
if ( ${zstarwinbad} == 1 ) then
  echo "#WARNING: Warning: zstarwin is ${zstarwin}, but reasonable would be ${zstarwintmp}" >> LOGS/${scriptname}.results
  echo "::WARNING: zstarwin is ${zstarwin}, but reasonable would be ${zstarwintmp}"
endif
#
echo "<<@progress: 1>>"
#
#################################################################################
${proc_2dx}/linblock "Verifying some parameters"
#################################################################################
#
if ( ${merge_alsoevenodd} == "y" && ( ! -e 2dx_merge_dirfile_even.dat || ! -e 2dx_merge_dirfile_odd.dat )) then
  echo ":: "
  echo ":: "
  ${proc_2dx}/linblock "ERROR: First define EVEN and ODD image tags, using the Project Tool"
  ${proc_2dx}/protest "This is the button top left in the GUI"
endif
#
if ( ${MergeHKMAX}x == 'x' ) then
  set MergeHKMAX = '20'
  ${proc_2dx}/linblock "ERROR: correcting MergeHKMAX to ${MergeHKMAX}"
  echo "set MergeHKMAX = ${MergeHKMAX}" >> LOGS/${scriptname}.results
endif
#
if ( `echo ${RESMIN} ${RESMAX} | awk '{ if ( $1 < $2 ) { s = 1 } else { s = 0 }} END { print s }'` == 1 ) then
  set oldval = ${RESMIN}
  set RESMIN = ${RESMAX}
  set RESMAX = ${oldval}
  ${proc_2dx}/linblock "ERROR: exchanging RESMIN and RESMAX, to RESMIN=${RESMIN}, and RESMAX=${RESMAX}"
  echo "set RESMIN = ${RESMIN}" >> LOGS/${scriptname}.results
  echo "set RESMAX = ${RESMAX}" >> LOGS/${scriptname}.results
endif
#
if ( ${ILIST}x == 'x' ) then
  set ILIST = "n"
  echo "set ILIST = ${ILIST}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "WARNING: ILIST corrected to ${ILIST}"
endif
#
if ( ${ILIST} == "n" ) then
  set ILIST_VAL = 0
else
  set ILIST_VAL = 1
endif
#
# This translates the list of directories to work on into one single long line:
cat 2dx_merge_dirfile.dat | tr "\n" " " > SCRATCH/2dx_merge_dirfile_oneline.dat
set dirlist = "`cat SCRATCH/2dx_merge_dirfile_oneline.dat`"
#
# This memorizes the current merge directory under the variable "olddir":
set olddir = $PWD
#
echo "The current working directory is" ${olddir}
#
echo "<<@progress: 5>>"
#
#############################################################################
${proc_2dx}/linblock "Sourcing sym2spsgrp_sub.com"
#############################################################################
#
source ${proc_2dx}/2dx_sym2spcgrp_sub.com
#
echo SYM = ${SYM}
echo spcgrp = ${spcgrp}
echo CCP4_SYM = ${CCP4_SYM}
#
############################################################################# 
${proc_2dx}/lin "2dx_merge_makedirs - to create all required subdirectories"
#############################################################################
#
source ${proc_2dx}/2dx_merge_makedirs
#
echo "<<@progress: 10>"
#
set create_merged_dataset = "y"
#
set NPRG = 1
if ( ${merge_reference} == '0' ) then
  # Use merge.aph files
  if ( ! -e APH/merge.aph ) then
    ${proc_2dx}/linblock "ERROR: APH/merge.aph not found. No reference for refinement available."
    ${proc_2dx}/protest "Run Merging script first."
  endif
  if ( ${reftiltgeo} == "y" ) then
    # tilt geometry refinement only possible with refined lattice lines in MTZ data, not with APH data.
    ${proc_2dx}/linblock "ERROR: tilt geometry refinement only possible with refined lattice lines in MTZ file."
    ${proc_2dx}/linblock "ERROR: Conflict between parameter USE MERGED DATA... and REFINE TILTGEOMETRY..."
    ${proc_2dx}/protest "If you want to refine the tilt geometry, switch to MTZ data first."
  endif
  set reference_file="APH/merge.aph"
endif 
if ( ${merge_reference} == '1' ) then
  # Use interpolated lattice lines merge3Dref_MRClefthanded.mtz
  if ( ! -e merge3Dref_MRClefthanded.mtz ) then
    ${proc_2dx}/linblock "ERROR: merge3Dref_MRClefthanded.mtz not found. No reference for refinement available."
    ${proc_2dx}/protest "Run Merging script first."
  endif
  set reference_file="merge3Dref_MRClefthanded.mtz"
  set NPRG = 3
endif

set itnum = 0

while ( ${itnum} < ${merge_refine_iterations} ) 
  #
  set itnum_p1 = ${itnum}
  @ itnum_p1 += 1
  echo ":: "
  echo "::         ITERATION ${itnum_p1} of ${merge_refine_iterations}"
  echo ":: "
  #
  set prog_num = `echo ${itnum} ${merge_refine_iterations} | awk '{ s = 10 + int( 80 * $1 / $2 ) } END { print s }'` 
  echo "<<@progress: ${prog_num}>>"
  #
  if ( ${merge_alsoevenodd} == "y" ) then
    echo ":: "
    #############################################################################
    ${proc_2dx}/linblock "Switching to EVEN IMAGES series"
    #############################################################################
    set dirfile = "2dx_merge_dirfile_even.dat"
    if ( ${merge_reference} == '0' ) then
      set reference_file="APH/merge_even.aph"
    endif 
    if ( ${merge_reference} == '1' ) then
      set reference_file="merge3Dref_MRClefthanded_even.mtz"
    endif  
    #
    ${proc_2dx}/linblock "Sourcing 2dx_origtilt_refine.com"
    source ${proc_2dx}/2dx_origtilt_refine.com
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_merge_generate3D.com"
    source ${proc_2dx}/2dx_merge_generate3D.com
    #
    cp APH/merge.aph APH/merge_even.aph
    cp APH/latfitted.hkl APH/latfitted_even.hkl
    cp merge3Dref_MRClefthanded.mtz merge3Dref_MRClefthanded_even.mtz
    #
    echo ":: "
    #############################################################################
    ${proc_2dx}/linblock "Switching to ODD IMAGES series"
    #############################################################################
    set dirfile = "2dx_merge_dirfile_odd.dat"
    if ( ${merge_reference} == '0' ) then
      set reference_file="APH/merge_odd.aph"
    endif 
    if ( ${merge_reference} == '1' ) then
      set reference_file="merge3Dref_MRClefthanded_odd.mtz"
    endif  
    #
    ${proc_2dx}/linblock "Sourcing 2dx_origtilt_refine.com"
    source ${proc_2dx}/2dx_origtilt_refine.com
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_merge_generate3D.com"
    source ${proc_2dx}/2dx_merge_generate3D.com
    #
    cp APH/merge.aph APH/merge_odd.aph
    cp APH/latfitted.hkl APH/latfitted_odd.hkl
    cp merge3Dref_MRClefthanded.mtz merge3Dref_MRClefthanded_odd.mtz
    #
    echo ":: "
    #############################################################################
    ${proc_2dx}/linblock "Switching to ALL IMAGES series"
    #############################################################################
    set dirfile = "2dx_merge_dirfile.dat"
  else
    set dirfile = "2dx_merge_dirfile.dat"
    if ( ${merge_reference} == '0' ) then
      set reference_file="APH/merge.aph"
    endif 
    if ( ${merge_reference} == '1' ) then
      set reference_file="merge3Dref_MRClefthanded.mtz"
    endif
    #
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_origtilt_refine.com"
    #############################################################################  
    source ${proc_2dx}/2dx_origtilt_refine.com
    #
  endif
  #
  if ( ${merge_alsoevenodd} == "y" ) then
    echo ":: "
    #############################################################################
    ${proc_2dx}/linblock "Switching to EVEN IMAGES series"
    #############################################################################
    set dirfile = "2dx_merge_dirfile_even.dat"
    #
    ${proc_2dx}/linblock "Sourcing 2dx_origtilt_merge.com"
    source ${proc_2dx}/2dx_origtilt_merge.com
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_merge_generate3D.com"
    source ${proc_2dx}/2dx_merge_generate3D.com
    #
    cp APH/merge.aph APH/merge_even.aph
    cp APH/latfitted.hkl APH/latfitted_even.hkl
    cp merge3Dref_MRClefthanded.mtz merge3Dref_MRClefthanded_even.mtz
    #
    echo ":: "
    #############################################################################
    ${proc_2dx}/linblock "Switching to ODD IMAGES series"
    #############################################################################
    set dirfile = "2dx_merge_dirfile_odd.dat"
    #
    ${proc_2dx}/linblock "Sourcing 2dx_origtilt_merge.com"
    source ${proc_2dx}/2dx_origtilt_merge.com
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_merge_generate3D.com"
    source ${proc_2dx}/2dx_merge_generate3D.com
    #
    cp APH/merge.aph APH/merge_odd.aph
    cp APH/latfitted.hkl APH/latfitted_odd.hkl
    cp merge3Dref_MRClefthanded.mtz merge3Dref_MRClefthanded_odd.mtz
    #
    echo ":: "
    #############################################################################
    ${proc_2dx}/linblock "Switching to ALL IMAGES series"
    #############################################################################
  else
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_origtilt_merge.com"
    #############################################################################
    #
    source ${proc_2dx}/2dx_origtilt_merge.com
    #
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_plttilt.com"
    #############################################################################
    #
    source ${proc_2dx}/2dx_plttilt.com
    #
    #############################################################################
    ${proc_2dx}/linblock "Sourcing 2dx_merge_generate3D.com"
    ############################################################################# 
    #
    source ${proc_2dx}/2dx_merge_generate3D.com
    #
  endif
  #
  # Issue evaluate here to reload the Ph-Ori-Change column
  echo "<<@evaluate>>"
  #
  @ itnum += 1
end
#
#
${proc_2dx}/linblock "Now, press the REFRESH DISPLAY VIEW button in top menu to update the project library."
#############################################################################
${proc_2dx}/linblock "2dx_refine normal end."
#############################################################################
#
echo "<<@progress: 100>>"
#
exit
# for GUI:

